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Abstract 



Methods of Computational Fluid Dynamics are applied 
to simulate pulsatile blood flow in human vessels and in 
the aortic arch. The non-Newtonian behaviour of the hu- 
man blood is investigated in simple vessels of actual size. 
A detailed time -dependent mathematical convergence test 
has been carried out. The realistic pulsatile flow is used 
in all simulations. Results of computer simulations of the 
blood flow in vessels of two different geometries are pre- 
sented. For pressure, strain rate and velocity component 
distributions we found signiflcant disagreements between 
our results obtained with realistic non-Newtonian treatment 
of human blood and widely used method in literature: a 
simple Newtonian approximation. A signiflcant increase of 
the strain rate and, as a result, wall sear stress distribution, 
is found in the region of the aortic arch. We consider this 
result as theoretical evidence that supports existing clini- 
cal observations and those models not using non-Newtonian 
treatment underestimate the risk of disruption to the human 
vascular system. 

Keywords: Fluid dynamics, Navier-Stokes equation, 
human blood flow, non-Newtonian viscosity, human ves- 
sels, aortic arch, pressure and wall shear stress distributions. 



1. Introduction 

Human blood is a liquid with variable density and vis- 
cosity. The movement of the blood inside vessels and ar- 
teries can be described by fundamental laws of physics, i.e. 
equations of fluid dynamics. The scientific literature now 
contains many citations where researchers have used com- 
puter simulations of blood flows in various size vessels and 
arteries at different spatial geometries, see for example [1- 
10]. Whereas experimental investigations of vascular dy- 
namics and flow are complicated or simply impossible to 
carry out due to small sizes of vessels in living systems. 
Therefore, theoretical-mathematical models and computer 
simulations are very useful for studying blood flows. 

Atherosclerosis occurs at specific arterial sites. This 
phenomenon is related to hemodynamics and to wall shear 
stress (WSS) distributions. WSS is the tangential drag force 
produced by moving blood. It is a mathematical function of 
the velocity gradient of blood near the endothelial surface: 
T„ = ^ [dU{t, y, Rv)/dy]y^Q here ^ is the dynamic vis- 
cosity, t is current time, U (t, y, Ry) is the flow velocity par- 
allel to the wall, y is the distance to the wall of the vessel, 
and Ry is its radius. It was shown, that the magnitude of 
WSS is directly proportional to blood flow and blood vis- 
cosity and inversely proportional to the cube of the radius 
of vessels, in other words a small change of the radius of a 



vessel will have a large effect on WSS. 

Arterial wall remodeling is regulated by WSS. In re- 
sponse to high shear stress arteries enlarge. Consequently, 
the atherosclerotic plaques localize preferentially in regions 
of low shear stresses, but not in regions of higher shear 
stresses. Furthermore, decreased shear stress induces inti- 
mal thickening in vessels which have adapted to high flow. 
Also, final vascular events that induce fatal outcomes, such 
as acute coronary syndrome, are triggered by the sudden 
mechanical disruption of an arterial wall. Thus, we can con- 
clude, that the final consequences of tragic fatal vascular 
diseases are strongly connected to mechanical events that 
occur on the vascular wall, and these, in turn, which are 
Ukely to be heavily influenced by alterations in blood flow 
and the characteristics of the blood itself. 

In order to predict, diagnose, and prevent fatal outcomes 
in these vascular diseases, fluid + solid mechanical inter- 
actions between the human blood and the vascular wall are 
attractive and necessary targets for analysis. However, it 
is very difficult to make measurements of detailed mechan- 
ical properties in living systems. In turn, theoretical bio- 
mathematics provides a series of logic tools that can over- 
come limitations of direct measurements in highly labile 
Uving systems and provide a framework for testing vari- 
ables, and is more rapid, efficient, and possibly more pre- 
dictive than repeated experimentation in animal modeUng 
systems 

In this work we carry out real-time full-dimensional 
computer simulations of pulsatile blood flows in actual size 
vessels and in the aortic arch. We take into account dif- 
ferent physical effects on blood flow and examine the non- 
Newtonian nature of human blood as a fluid. 

Computer simulation is well-suited to those cases in 
which it is difficult to carry out reliable experiments due 
to the very small size of some vessels, such as, for example, 
coronary arteries. Also, bio-mathematical computer simu- 
lations may be especially useful in situations when an in- 
travascular stent is implanted inside a vessel, aortic arch, 
aortas etcetera. Stents may have compUcated shapes and 
they are very small micro-devices (typically only several 
millimeters in diameter, when fuUy deployed). It is use- 
ful to know pressure, strain rate distributions and the profile 
of blood flow after coming through a stented vessel. 

The next section presents the mathematical and numer- 
ical methods used in this work. Section 3 presents our re- 
sults. The CGS unit system is used in all simulations as 
well as for presentation of our results. 

2. Method 

In this work we consider the human blood as an incom- 
pressible fluid, and flows in vessels and in the aortic arch 
are assumed to be laminar. One of the goals of this work 



is to investigate the non-Newtonian behaviour of the human 
blood. To attain these ends we carry out simulations for 
same systems, but with different models of the blood. Then 
we compare the results. 

For many years, investigation of hemorheology has been 
of great interest in the field of biomedical engineering. 
Researchers have investigated correlations for example, of 
stroke, arterial diseases, hypertension, and whole blood vis- 
cosity. Blood consists of plasma and particles, including 
red blood cells, leukocytes, platelets and macromolecular 
protein aggregates. The viscosity of blood depends on the 
viscosity of plasma, in combination with the hematocrit (a 
measure of the particulate component of blood). The nor- 
mal hematocrit of human blood ranges between 35% - 45%. 
Higher hematocrit implies higher viscosity. The relation be- 
tween hematocrit and viscosity is very complex and in the 
scientific hterature, many mathematical fitting formulas are 
available for assessing this relationship. 

Next, the viscosity of blood determines its velocity. That 
is, when velocity or shear rate increases viscosity decreases. 
Also, the viscosity/velocity depends on the size of the blood 
vessel. This is called the Fahraeus-Lindqvist effect, that 
is in small diameter blood vessels, and at higher veloci- 
ties, blood viscosity decreases. Viscosity of human blood 
strongly depends on its temperature. 

To carry out our simulations we used a commercial pro- 
gram FLOW3D from Flow Science Inc., Santa Fe, New 
Mexico, USA. FLOW3D is a general purpose CFD pack- 
age. It appUes speciaUy developed numerical techniques to 
solve the equations of motion of fluid. The methods imple- 
mented in the program allow us to obtain transient, 3D solu- 
tions for multi-scale and multi-physics flow problems. One 
of the most attarctive sides of FLOW3D is very well de- 
veloped numerical techniques to solve non-linear fluid dy- 
namic equations. To our best knowledge our work is the first 
time application of the FLOW3D program to blood flows in 
human vessels. 

When the turbulence option is used, the viscosity is 
a sum of the molecular and turbulent values. For non- 
Newtonian fluids the viscosity can be a function of the strain 
rate and/or temperature. A general expression based on the 
Carreau model is used in FLOW-3D for the strain rate de- 
pendent viscosity: 

_ MO-E-T — Moo 

" ^ Aoo + [Ao + (Aiii;T)2ei,ey](i-")/2 
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where = l/2{dui/ dxj + duj/dxi) is the fluid strain 
rate in Cartesian tensor notations, ^oo, MOj Aq, Ai, A2 and n 
are constants. Also, Et = exp[a{T* / {T - b) - C)], where 
T*,a,b, and c are also parameters of the temperature de- 
pendence, and T is fluid temperature. This basic formula is 



used in our simulations for blood flow in vessels and in the 
aortic arch. 

Next, fluid dynamics is described with 2-nd order non- 
linear, transient differential equations. The governing equa- 
tions consist of the continuity equation and the Navier- 
Stokes equations. The general mass continuity equation, 
which is solved within the FLOW3D program has the fol- 
lowing general form: 

= R^f + Rsor, (2) 

X 

where Vp is the fractional volume open to flow, R and xi 
are coefficients which value dependent on the coordinate 
system: {x,y,z) or {r,d,z), p is the fluid density, 
is a turbulent diffusion term, and Rsor is a mass source, 
(m, V, w) are the velocity components in coordinate direc- 
tions {x,y,z) respectively. For example, when Cartesian 
coordinates are used, R — 1 and ^ = 0, see FLOW3D man- 
ual Q. Finally, is the fractional area open to flow in the 
X direction, analogously for Ay and A^. 
The turbulent diffusion term is 

d . .dp. d . . r^dp. 

+ (3) 

oz oz X 

where the coefficient Vp = Cpp,/ p, p. is dynamic viscosity 
and Cp is a constant. The Rsor term is a density source term 
that can be used to model mass injections through porous 
obstacle surfaces. 

It is well known, that compressible flow problems re- 
quire solution of the full density transport equation. In this 
work we treat blood as an incompressible fluid. For incom- 
pressible fluids p = constant and the Eqn. (2) becomes the 
following: 

R. 



(4) 



The equations of motion for the fluid velocity compo- 
nents (m, u, w) in the 3-coordinate system are the Navier- 
Stokes equations with specific additional terms included in 
the FLOW3D program: 

du ^ I j\ A '^^ + A '^^ ^^^^^ 
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where, {Gx, Gy, Gz) are so called body accelerations ||9l, 
{fx,fy,fz) are viscous accelerations, {bx,by,bz) are flow 
losses in porous media or across porous baffle plates, and 
the final term accounts for the injection of mass at a source 
represented by a geometry component. As we mentioned 
above, FLOW3D is a general purpose fluid dynamics pro- 
gram, which includes many specific situations. However, in 
this short description we try to make a valuable impression 
about FLOW3D and give here the general form of all equa- 
tions. Next, the term — (u^,, w^,, w^y) in Eqn. (5) is 
the velocity of the source component, which will generally 
be non-zero for a mass source at a General Moving Object 
(GMO) Jo). The term Us = {us, Vs,Ws) is the velocity of 
the fluid at the surface of the source relative to the source 
itself. It is computed in each control volume as 



Us 



1 d{Qn) 
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(8) 



where dQ is the mass flow rate, ps fluid source density, dA 
the area of the source surface in the cell and n the outward 
normal to the surface. 

The source is of the stagnation pressure type when in 
Eqs. dSllTl l S = 0.0. Next, 6 = 1.0 corresponds to the 
source of the static pressure type. 

It is assumed, that at a stagnation pressure source fluid 
enters the domain at zero velocity. As a result, pressure 
should be considered at the source to move the fluid away 
from the source. For example, such sources are designed to 
model fluid emerging at the end of a rocket or the simple de- 
flating process of a balloon. In general, stagnation pressure 
sources apply to cases when the momentum of the emerg- 
ing fluid is created inside the source component, like in a 
rocket engine. At a static pressure source the fluid velocity 
is computed from the mass flow rate and the surface area of 
the source. In this case, no extra pressure is required to pro- 
pel the fluid away from the source. An example of such a 
source is fluid emerging from a long straight pipe. Note that 
in this case the fluid momentum is created far from where 
the source is located. 

For a variable dynamic viscosity p, the viscous acceler- 
ations are 
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Figure 1. Velocity waveform at the vessel in- 
let. Results taken from Fig. 2 of work t6j. 
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Wy and wl are wall shear 



In Eqs. (fTZt-dTTIi the terms w: 
stresses. If these terms are equal to zero, there is no wall 
shear stress. This is because the remaining terms contain 
the fractional flow areas {A^, Ay, A^) which vanish at the 
walls. The wall stresses are modeled by assuming a zero 
tangential velocity on the portion of any area closed to flow. 
Mesh and moving obstacle boundaries are an exception be- 
cause they can be assigned non-zero tangential velocities. 
In this case the allowed boundary motion corresponds to a 
rigid body translation of the boundary parallel to its surface. 
For turbulent flows, a law-of-the-wall velocity profile is as- 
sumed near the wall, which modifies the wall shear stress 
magnitude. 

As we already mentioned, in all simulations we consider 
the blood flow as pulsatile flow. The final result for the 
inflow waveform has been taken from work |6l. The wave- 
form is shown in Fig. 1. These velocity values are used as 
time-dependent inflow initial boundary conditions. These 
numbers are included directly in the FLOW3D program. 

The equations of fluid dynamics should be solved to- 
gether with specific boundary conditions. The numerical 



model starts with a computational mesh, or grid. It consists 
of a number of interconnected elements, or 3D-cells. These 
3D-cells subdivide the physical space into small volumes 
with several nodes associated with each such volume. The 
nodes are used to store values of the unknown parameters, 
such as pressure, strain rate, temperature, velocity compo- 
nents etcetera. This procedure provides values for defining 
the flow parameters at discrete locations and to set up spe- 
cific boundary conditions. Finally, one can start developing 
effective numerical approximations for the solution of the 
fluid dynamics equations. 

New pressure-velocity solvers have been implemented 
in FLOW-3D. We used the so called GMERS method. GM- 
RES stands for the generalized minimum residual method. 
In addition to the GMRES solver, a new optional algorithm, 
the generalized conjugate gradient (GCG) algorithm, has 
also been implemented for solving viscous terms in the new 
GMRES solver. This new solver is a highly accurate and 
efficient method for a wide range of problems. It possesses 
good convergence, symmetry and speed properties; how- 
ever, it does use more memory than the SOR or SADI meth- 
ods. The GMRES solver does not use any over- or under- 
relaxation ||9|. 

3. Results 

This section represents the results of our simulations for 
blood flows in a simple vessel and in the human aortic arch. 
The first geometry was chosen for preliminary test calcula- 
tions and testing the FLOW3D program. One of the most 
important preliminary testing tasks is the check numerical 
convergence. This test has been successful and our results 
will be shown below in this paper. To our best knowledge 



the current work is a first attempt to apply FLOW3D to hu- 
man blood flows in vessels and in aortic arch. 

First, we present results for a simpler geometry vessel in 
the shape of a tube. However, the human blood is treated 
as real and a non-Newtonian liquid. The necessary data for 
viscosity of the blood we found from previous laboratory 
and clinical measurements lHJO. We take into account the 
real pulsatile flow, which is shown in Fig. 1. The data for 
Fig. 1 have also been obtained in clinical measurements ||6l. 

After such preliminary simulations we switch to a more 
complicated spatial configuration. In this work it is the aor- 
tic arch. It is axiomatic that real people may have different 
size aortic arches with slightly different shapes. However, 
we carried out simulations for an average size and shape 
aortic arch. 

The main goal of this work is to treat the above men- 
tioned systems realistically, reveal the physics of the blood 
flow dynamics, and to obtain reliable results for pressure, 
dynamic viscosity, velocity profiles and strain rate distribu- 
tions. Also, we tested, the widely cited in literature, New- 
tonian and non-Newtonian models of the human blood. 

3.1 Blood flow in vessel 

As we mentioned above in this work we adopted the 
shape of a straight vessel as a tube. The sizes of the tube 
are: L = 8 cm in length and R = 0.34 cm in the inner 
radius. The thickness of the vessel wall is s = 0.03 cm. We 
have chosen 5.5 cycles of the blood pulse. 

Consider in more detail the expression ([T]l. In these cal- 
culations we follow the works fH O, where the Carreau 
model of the human blood has also been used. In con- 
sistence with lU |3| we choose the following coefficients: 
A2 = Aoo = 0, a = and Et ~ 1, that is we don't take 
into account the temperature dependence of the viscosity. 
Next: Ao = 1, Ai = 3.313 sec, ^00 = 0.0345 P, /io = 0.56 
P, and n = 0.3568. 

In our calculations we applied a cylindrical coordinate 
system with the axis OZ directed over the tube axis. Dif- 
ferent numbers of cells have been used to discretize the 
empty space inside the tube. In the open space (inner part of 
the tube) the fluid dynamics equations have been solved to- 
gether with appropriate mathematical boundary conditions. 
The convergence was achieved when we used 52,800 cells, 
that is we used 100 points over OZ, 22 points over the ra- 
dius of the inside space R — 0.34 cm, and 24 points over 
azimuthal angle <i> from to 2tt. 

Time-dependent results for pressure, strain rate and ve- 
locity component W are presented in Fig. 2. We chose 
three precise geometrical points to compare results - 1st: in 
the inlet of the tube, when Z = 0; 2nd: in the middle point 
Z — 4.0 cm, and 3rd: in outlet point z — 8.0 cm. The data 
for Fig. 2 was obtained with two different models of human 



blood. The bold lines are results with the non-Newtonian 
viscosity ([T]) and the dashed lines are results with the New- 
tonian model when the viscosity /i has a constant value and 
equals to 0.0345 P. As one can see the results are different 
for strain rate distributions and very different for pressure 
distributions. These results clearly indicate that probably in 
most cases of computer simulations of human blood flows 
only the non-Newtonian model should be used. 

In Fig. 3 we show our time-dependent convergence 
results. These data are obtained with the realistic non- 
Newtonian model of the blood. Again, we have chosen 
three geometrical points to compare results: 1- in the in- 
let of the tube, when Z = 0; 2- in the middle point Z = 4.0 
cm, and 3- outlet point z — 8.0 cm. 

Three upper rows represent time-dependent plots for 
pressure, dynamics viscosity and strain rate distributed over 
the simulation time, which is 5.5 sec. However, below in 
three lower rows results for the kinematic characteristics of 
the blood flow are also shown, which are three components 
of the blood velocity: U, V, W. The results are obtained 
with three sets of computation cell distributions, where we 
used 36,000 cubic cefls, 47,000 and 52,000 cubic cefls. As 
we see from Fig. 3 it is harder to obtain convergence for 
the pressure distributions and easier for other shown param- 
eters. 

3.2 Blood flow in aortic arch 

The aortic arch is represented as a curved tube H. In our 
simulations the outer radius of the tube is 2.6 cm. A straight 
vessel (tube) is also merged to the arch. The length of the 
straight tube is about 4 cm. Again, the thickness of the wall 
is 0.03 cm, and the inner radius of the tube is r = 0.34 cm. 
The geometry is shown in Figs. 4 and 5. This configuration 
closely models and represents the real aortic arch. One of 
the goals of these simulations is to reveal the physics of the 
blood flow dynamics in the arch. 

Now we use the Cartesian coordinate system. Here 
we also carried out a convergence test. To better repre- 
sent the shape of the arch we applied five Cartesian sub- 
coordinate systems in our FLOW3D simulations. After the 
discretization the total number of all cubic cells reached 
about 400,000. It is important to mention here, that we 
again obtained a full numerical convergence. 

As we mentioned the goal of these simulations is to com- 
pute pressure, velocity and strain rate distributions in the 
arch, while the human blood is treated as a non-Newtonian 
liquid and the realistic pulsatile blood flow is used as it is 
shown in Fig. 1. In Figs. 4 and 5 we show the results 
for strain rate distributions inside the arch for four specific 
time moments. At the most left point, which is inlet, we 
specify the pulsatile velocity source as the initial condition, 
that is the data from Fig. 1 are used. From the general 
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Figure 2. Time-dependent pressure, strain 
rate and velocity component W - over OZ 
axis. Results for straight vessel for 3 spatial 
points: Z=0, 4, 8 cm. Bold lines: calculations 
with realistic non-Newtonian viscosity of hu- 
man blood. Dashed lines: simple Newtonian 
approximation. 



theory of fluid mechanics fTOl it is possible to determine 
together with viscosity and spatial geometry, the dynamics 
of the blood according to the Navier-Stokes equation and 
its boundary conditions. Small vectors indicate the blood 
velocity. As can be seen blood flows from left to right in 
direction. However, because of pulsatility blood flows in 
the opposite direction too. It is seen in Fig. 4 - lower right 
graph #43. It is in good agreement with the general physi- 
cal intuition, and it additionally shows correctness of these 
simulations. 

The values of the strain rate are also shown. These values 
are strongly oscillating. From the plots one can conclude 
that in the region of the arch the strain rate values are be- 
coming much larger than in the region of the straight vessel. 
This result represents clear evidence that in this part of the 
human vascular system atherosclerotic plaques should lo- 
calize less than in the straight vessels. However, the higher 
wall shear stress values in the aortic arch could be the rea- 
son for sudden mechanical disruption of the arterial wall in 
this part of the human vascular system. These results are 
consistent with laboratory and clinical observations. 



In conclusion, we would like to point out here, that the 
developments in this work can be directly applied to even 
more interesting and important situation such as, when a 
stent is implanted inside a vessel |7] [H. In that case, 
for example, it is important to determine blood flow dis- 
turbance, the pressure distribution, strain rate and etcetera. 
This work is in progress in our group. 
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Figure 3. Time-dependent convergence test: results for pressure, dynamic viscosity, strain rate 
and velocity components: U over OX, V over OY, and W over OZ. 



Figure 4. Blood flow in the aortic arch for two consecu- 
tive moments of the discretized time. A strong pulsatility 
of the strain rate values is seen: Upper plot shows results 
for t = 4.218 sec, where strain rate ranges from e = 1.1 
1/sec to e = 24.6 1/sec. Lower plot shows results for t = 
4.329 sec, where strain rate ranges from e = 5. 1/sec to 
e=234. 1/sec. The maximum values of the strain rate are 
localized in the region inside the arch. Blood flows from 
right to left in both pictures. 



Figure 5. Blood flow in the aortic arch for two consecu- 
tive moments of discretized time. Upper plot shows re- 
sults for t = 4.551 sec, where strain rate ranges from e = 
0.4 1/sec to e = 22.6 1/sec. Lower plot shows results for t 
= 4.662 sec, where strain rate ranges from e = 0.1 1/sec 
to e = 27.6 1/sec. The maximum values of the strain rate 
are localized again in the region inside the arch. In the 
upper plot blood flows from left to right, however in the 
lower plot: from right to left. 



